function bspline2

%  Plots the error and cpu time using FDM and B-splines 
%  for the BVP
%       y'' + p(x)y' + q(x)y= f(x)   for xL < x < xr  '
%  where 
%      y(xl) = yL  and y(xR) = yR

% clear all previous variables and plots
clear *
clf

% set boundary conditions
	xL=0; yL=0;
	xR=1; yR=exp(-2);

% determine error and cpu time as number of points increases
ii=0;
for i=0:2
	n=4*10^i+1;
	ii=ii+1; points(ii)=n;
	x=linspace(xL,xR,n);
	tic
	y=FDM(x,yL,yR,n);
	time_fdm(ii)=toc;
	tic
	yy=BSM(x,yL,yR,n);
	time_bsm(ii)=toc;
	for k=1:n
		exact(k)=x(k)*exp(-2*x(k));
	end;
	errorM(ii)=norm(exact-y,inf);
	errorMM(ii)=norm(exact-yy,inf);
end;

% plot error values
subplot(2,1,1)
loglog(points,errorM,'--rs','LineWidth',1,'MarkerSize',7)
hold on
loglog(points,errorMM,'--bo','LineWidth',1,'MarkerSize',7)
legend(' FDM',' B-spline',1)
grid on
set(gca,'MinorGridLineStyle','none')

% define title and axes used in plot
xlabel('Number of Grid Points','FontSize',14,'FontWeight','bold')
ylabel('Error','FontSize',14,'FontWeight','bold')
title('B-Splines vs FDM','FontSize',14,'FontWeight','bold')

% have MATLAB use certain plot options (all are optional)
box on
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 
% Set legend font to 14/bold                            		
set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold'); 
axis([1 1e3 1e-7 1e-1]);
set(gca,'ytick',[1e-9  1e-7 1e-5 1e-3 1e-1]);
hold off

% plot time values
subplot(2,1,2)
loglog(points,time_fdm,'--rs','LineWidth',1,'MarkerSize',7)
hold on
loglog(points,time_bsm,'--bo','LineWidth',1,'MarkerSize',7)
legend(' FDM',' B-spline',3)
grid on
set(gca,'MinorGridLineStyle','none')

% define title and axes used in plot
xlabel('Number of Grid Points','FontSize',14,'FontWeight','bold')
ylabel('CPU Time','FontSize',14,'FontWeight','bold')

% have MATLAB use certain plot options (all are optional)
box on
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 
% Set legend font to 14/bold                            		
set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold'); 
axis([1 1e3 1e-4 1e-1]);
set(gca,'ytick',[1e-4  1e-3 1e-2 1e-1 1e0]);
hold off


function g=q(x)
g=-2;

function g=p(x)
g=1;

function g=f(x)
g=-3*exp(-2*x);


% finite difference method	
function y=FDM(x,yL,yR,n)
h=x(2)-x(1);
a=zeros(1,n-2); b=zeros(1,n-2); c=zeros(1,n-2);
for i=1:n-2
	a(i)=-2+h*h*q(x(i+1));
	b(i)=1-0.5*h*p(x(i+1));
	c(i)=1+0.5*h*p(x(i+1));
	z(i)=h*h*f(x(i+1));
end;
z(1)=z(1)-yL*b(1);
z(n-2)=z(n-2)-yR*c(n-2);
y=tri(a,b,c,z);
y=[yL, y, yR];


% B-spline method
function y=BSM(x,yL,yR,N)
h=x(2)-x(1);
a=zeros(1,N); b=zeros(1,N); c=zeros(1,N); z=zeros(1,N);
hh=h*h;
for i=1:N
	a(i)=4*(-3+hh*q(x(i)));
	b(i)=6-3*h*p(x(i))+hh*q(x(i));
	c(i)=6+3*h*p(x(i))+hh*q(x(i));
	z(i)=6*hh*f(x(i));
end;
z(1)=z(1)-6*yL*b(1); a(1)=a(1)-4*b(1);  c(1)=c(1)-b(1);
z(N)=z(N)-6*yR*c(N); a(N)=a(N)-4*c(N);  b(N)=b(N)-c(N);
w=tri(a,b,c,z); 
wL=6*yL-4*w(1)-w(2); wR=6*yR-w(N-1)-4*w(N);
ww=[wL, w, wR];
y(1)=yL;
for ix=2:N-1
	y(ix)=(ww(ix)+4*ww(ix+1)+ww(ix+2))/6;
end;
y(N)=yR;


function y=bspline(x)
% Calculate the value of a cubic B-spline at point x
x=abs(x) ;
if x > 2,
  y=0 ;
else
  if x > 1,
    y=(2-x)^3/6 ;
  else
    y=2/3-x^2*(1-x/2) ;
  end ;
end ;


% tridiagonal solver
function y = tri( a, b, c, f )
N = length(f);
v = zeros(1,N);   
y = v;
w = a(1);
y(1) = f(1)/w;
for i=2:N
    v(i-1) = c(i-1)/w;
    w = a(i) - b(i)*v(i-1);
    y(i) = ( f(i) - b(i)*y(i-1) )/w;
end
for j=N-1:-1:1
   y(j) = y(j) - v(j)*y(j+1);
end